function [dof_boundary, ub_1, ub_2] = dirichlet_bb(dof_map, V, T, TE, ET, FE_Order, bdr, g, varargin)
% function [dof_boundary, ub_1, ub_2]  = dirichlet_bb(dof_map, V, T, TE, ET, FE_Order, bdr, g, varargin)
%
% Input: 
%   dof_map       ----- the dof mapping for each elements
%   (V, T, TE,ET) ----- the triangular mesh informations
%   FE_Order     ----- the order of polynomail spaces in BB Form
%   bdr               ----- the Dirichlet boundary edge list
%   g                  ----- the boundary function
%
% Output:
%   boundary_dof    ----- edge-wise boundary dof index ((FE_Order + 1) x n_bdr)
%   ub        ----- the value of boundary dof, in two dimension
%
%  Dr. Xian-Liang Hu
%  Aug 2012
%

% some constant
cr = cr_pattern(FE_Order);
n_bdr = length(bdr);
n_dof_per_edge = (FE_Order + 1);
dof_boundary = zeros(n_dof_per_edge, n_bdr);

v_start = zeros(n_bdr, 1);
v_end  =  zeros(n_bdr, 1);


for k = 1:n_bdr
    eg = bdr(k); tri = ET(eg,1);
    eg_local_idx = find(TE(tri,:)==eg);
    
    v_start_idx= mod(eg_local_idx, 3) + 1;
    v_end_idx = mod(v_start_idx, 3) + 1;
    v_start(k) = T(tri, v_start_idx);
    v_end(k)  =  T(tri, v_end_idx);
    
    if (eg_local_idx ~= 2)
        dof_local = cr_indices(0, FE_Order, eg_local_idx, cr);
    else
        dof_local = cr_indices(0, FE_Order, -eg_local_idx, cr);
    end
    dof_boundary(:, k) = dof_map(dof_local, tri);
end

px = (V(v_end, 1)*(0:FE_Order) + V(v_start, 1)*(FE_Order:-1:0))/FE_Order;
py = (V(v_end, 2)*(0:FE_Order) + V(v_start, 2)*(FE_Order:-1:0))/FE_Order;

[g_1, g_2] = feval(g, px', py', varargin{:});   

% additional treatment for BB polynomials
B_mat = vdm11_bb(FE_Order);  % the common interp matrix for all edges
ub_1 = B_mat\g_1;
ub_2 = B_mat\g_2;   % obtain the value of B-coefficient by solving interpolation e-b-e
end

%%%%%%%%%
%
%
function Mat = vdm11_bb(d)
I=d:-1:0;
J=0:d;
IM = diag(I)*ones(d+1,d+1);
JM = diag(J)*ones(d+1,d+1);
Mat = (IM/d).^(IM').*(JM/d).^(JM');
IF = gamma(I+1);
JF = gamma(J+1);
A = factorial(d)*ones(d+1,d+1)*diag(1./(IF.*JF));
Mat = A.*Mat;
end
